In this practical, we will explore our data in more detail using some basic statistical tests to identify where the variation in our data lies, and hence where the interesting ecological insights might lie. In this example dataset, we have d13C, d15N and C:N ratio from three fish species.
mydata <- read.csv("Practical02.csv", header = TRUE, stringsAsFactors = FALSE)
# force species to be a factor type variable
mydata$Species <- factor(mydata$Species)
# verify that it imported as expected
head(mydata)
NA
# code to do this using base R graphics
# plot(Liver.d13C ~ Liver.CN, data = mydata,
# subset = Species == "Arrowtooth.Flounder", main = "Arrowtooth.Flounder")
#
# # lowess smoother not working like this for the subset.
# # lines(lowess(mydata$Liver.CN, mydata$Liver.d13C), col = "red")
#
# plot(Liver.d13C ~ Liver.CN, data = mydata,
# subset = Species == "Pacific.cod", main = "Pacific.cod")
#
# plot(Liver.d13C ~ Liver.CN, data = mydata,
# subset = Species == "Pollock", main = "Pollock")
# try to do this in ggplot using facets
library(tidyverse)
liver.plot <- ggplot(data = mydata,
mapping = aes(x = Liver.CN,
y = Liver.d13C)) +
geom_point() +
geom_smooth(method = "loess", alpha = 0.5) +
facet_grid(Species~.) +
theme(text = element_text(size=20))
print(liver.plot)
`geom_smooth()` using formula 'y ~ x'
The liver in particular has a large quantity of lipids, and in this example, we apply a simple numerical correction based on Logan et al 2008. This correction is weighted by the C/N ratio, such that high amounts of C require a greater correction.
# constants used in the correction for liver tissue
a <- 6.059
b <- -22.270
c <- -1.397
# add columns for lipid corrected liver to the data.frame
mydata <- mydata %>%
mutate(Liver.d13C.crtd = Liver.d13C +
(a * Liver.CN + b) / (Liver.CN + c))
# mydata$Liver.d13C.crtd <- mydata$Liver.d13C +
# (a * mydata$Liver.CN + b) / (mydata$Liver.CN + c)
# check that it worked as expected
head(mydata)
# specify a color pallete with 3 colours for this example
palette(c("black", "red", "blue"))
# plot the correlation between the corrected and raw values, and
# indicate which observations have large C/N ratios by the
# size of the datapoints.
# plot(Liver.d13C.crtd ~ Liver.d13C, data = mydata,
# col = Species, cex = 5 * Liver.CN / max(Liver.CN), pch = 19, asp = 1)
#
# # add the 1:1 line
# abline(a = 0, b = 1, col = "grey", lty = 1, lwd = 2)
#
#
# # reset the colour palette to the default
# palette("default")
# now plot the corrected values on top of our raw data from before
# and add a vertical line at 3.6 which is "pure protein"
liver.plot.2 <- liver.plot +
geom_point(data = mydata, aes(x = Liver.CN,
y = Liver.d13C.crtd),
col = "red") +
geom_vline(xintercept = 3.6, col = "black")
print(liver.plot.2)
`geom_smooth()` using formula 'y ~ x'
As always, it is sensible to start by taking a quick look at the summary statistics describing the data, and how the data are encoded in R. The only bit of code I’ve changed here is to make the text of the correlation coefficient text a bit smaller when defining the function panel.cor().
# how are the data encoded in R?
str(mydata)
'data.frame': 113 obs. of 9 variables:
$ Species : Factor w/ 3 levels "Arrowtooth.Flounder",..: 1 1 1 1 1 1 1 1 1 1 ...
$ Length : int 320 170 180 360 160 500 420 350 370 480 ...
$ Muscle.d13C : num -19.1 -19.3 -19.3 -18.5 -19.5 ...
$ Muscle.d15N : num 15.3 13.3 12.8 15.3 13 ...
$ Muscle.CN : num 3.21 3.21 3.13 3.23 3.27 ...
$ Liver.d13C : num -22.5 -21.4 -22.2 -20.1 -22.6 ...
$ Liver.d15N : num 15.1 13 13.2 14.5 12.3 ...
$ Liver.CN : num 6.96 7.5 8.65 4.43 9.14 ...
$ Liver.d13C.crtd: num -18.9 -17.6 -18.1 -18.6 -18.4 ...
# a basic summary of each column in the data.frame
summary(mydata)
Species Length Muscle.d13C Muscle.d15N
Arrowtooth.Flounder:29 Min. : 90 Min. :-23.41 Min. :12.49
Pacific.cod :39 1st Qu.:270 1st Qu.:-19.75 1st Qu.:14.57
Pollock :45 Median :480 Median :-18.83 Median :16.03
Mean :450 Mean :-19.05 Mean :16.18
3rd Qu.:620 3rd Qu.:-18.11 3rd Qu.:17.79
Max. :820 Max. :-16.77 Max. :20.10
Muscle.CN Liver.d13C Liver.d15N Liver.CN
Min. : 3.105 Min. :-28.48 Min. : 9.896 Min. : 3.225
1st Qu.: 3.184 1st Qu.:-24.08 1st Qu.:13.336 1st Qu.: 6.395
Median : 3.235 Median :-22.71 Median :15.125 Median :10.189
Mean : 3.901 Mean :-22.62 Mean :15.032 Mean :15.793
3rd Qu.: 3.327 3rd Qu.:-21.38 3rd Qu.:16.750 3rd Qu.:23.063
Max. :21.149 Max. :-17.35 Max. :20.273 Max. :70.376
Liver.d13C.crtd
Min. :-22.85
1st Qu.:-19.03
Median :-18.41
Mean :-18.62
3rd Qu.:-17.92
Max. :-16.23
# how many obverations do we have per Taxon might also be of use
table(mydata$Species)
Arrowtooth.Flounder Pacific.cod Pollock
29 39 45
# Using datq from our first practical, generate all pair-wise scatterplots
# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
# source("my_functions.R")
# These functions are copied from the help file for pairs()
## put histograms on the diagonal
panel.hist <- function(x, ...)
{
usr <- par("usr"); on.exit(par(usr))
par(usr = c(usr[1:2], 0, 1.5) )
h <- hist(x, plot = FALSE)
breaks <- h$breaks; nB <- length(breaks)
y <- h$counts; y <- y/max(y)
rect(breaks[-nB], 0, breaks[-1], y, col = "cyan", ...)
}
panel.cor <- function(x, y, digits = 2, prefix = "", cex.cor, ...)
{
usr <- par("usr"); on.exit(par(usr))
par(usr = c(0, 1, 0, 1))
r <- cor(x, y)
txt <- format(c(r, 0.123456789), digits = digits)[1]
txt <- paste0(prefix, txt)
if(missing(cex.cor)) cex.cor <- 0.4/strwidth(txt)
# text(0.5, 0.5, txt, cex = cex.cor * r)
text(0.5, 0.5, txt, cex = 2)
}
# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
# this select function comes from the package dplyr and allows us to
# select a subset of columns to plot.
pairs(mydata,
diag.panel = panel.hist,
upper.panel = panel.smooth,
lower.panel = panel.cor)
NA
NA
NA
Using code from the first practical, we can use the package dplyr to generate summary statistics on the means and variation of each species in the dataset. I summarise the muscle and liver (using the lipid-corrected d13C) data separately only to keep the resulting table smaller and nicer to print.
# Summarise the data by group using the package dplyr (install if required)
# load the library
library(dplyr)
summarise.muscle <- mydata %>% group_by(Species) %>%
summarise(count = length(Species),
musC = mean(Muscle.d13C), musSdC = sd(Muscle.d13C),
musN = mean(Muscle.d15N), musSdN = sd(Muscle.d15N),
musCN = mean(Muscle.CN), musSdCN = sd(Muscle.CN))
# print a pretty table
knitr::kable(summarise.muscle, digits = 2)
| Species | count | musC | musSdC | musN | musSdN | musCN | musSdCN |
|---|---|---|---|---|---|---|---|
| Arrowtooth.Flounder | 29 | -20.08 | 1.30 | 14.64 | 1.2 | 4.30 | 2.21 |
| Pacific.cod | 39 | -18.07 | 1.07 | 17.98 | 1.0 | 3.78 | 2.41 |
| Pollock | 45 | -19.23 | 1.12 | 15.60 | 1.7 | 3.75 | 2.82 |
# print(summarise.muscle)
summarise.liver <- mydata %>% group_by(Species) %>%
summarise(count = length(Species),
livC.c = mean(Liver.d13C.crtd),
livSdC.c = sd(Liver.d13C.crtd),
livN = mean(Liver.d15N),
livSdN = sd(Liver.d15N),
livCN = mean(Liver.CN),
livSdCN = sd(Liver.CN) )
# print a pretty table
knitr::kable(summarise.liver, digits = 2)
| Species | count | livC.c | livSdC.c | livN | livSdN | livCN | livSdCN |
|---|---|---|---|---|---|---|---|
| Arrowtooth.Flounder | 29 | -18.66 | 0.70 | 14.21 | 1.34 | 8.38 | 3.79 |
| Pacific.cod | 39 | -18.02 | 1.18 | 16.77 | 1.48 | 12.65 | 7.97 |
| Pollock | 45 | -19.12 | 1.22 | 14.06 | 1.97 | 23.29 | 15.43 |
# print(summarise.liver)
The most suitable way to visualise the variation in the key tissues by species is probably to use simple boxplots. From here on in we will probably focus on the muscle data only and will put analysis of the liver samples on hold.
boxplot(Length ~ Species, data = mydata, col="grey",
xlab="Species", ylab="Length (cm)",
main = " variation in fish size")
boxplot(Muscle.d13C ~ Species, data = mydata, col="grey",
xlab="Species", ylab="Muscle d13C",
main = " variation in muscle d13C")
#now d15N
boxplot(Muscle.d15N ~ Species, data = mydata, col="grey",
xlab="Species", ylab="Muscle d15N",
main = " variation in muscle d15N")
#now d15N
boxplot(log(Muscle.CN) ~ Species, data = mydata, col="grey",
xlab="Species", ylab="log Muscle C/N",
main = " variation in muscle C/N")
NA
NA
An appropriate statistical test to determine whether the means of these boxplots are significantly different from one another is a simple ANOVA. You need to take a little care here, as the function you want when analysing data is called aov() in R (for Analysis of Variance) whereas there is also a function anova() but it is used to compare fitted linear model objects (a discussion point for another day, and one i wouldnt recommend in any case). I have log-transformed the C/N data owing to the heavy skew evident in the pairs plot above.
# length by species
length.aov <- aov(Length ~ Species, data = mydata)
summary(length.aov)
Df Sum Sq Mean Sq F value Pr(>F)
Species 2 572074 286037 7.872 0.000638 ***
Residuals 110 3996926 36336
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# muscle carbon by species
muscleC.aov <- aov(Muscle.d13C ~ Species, data = mydata)
summary(muscleC.aov)
Df Sum Sq Mean Sq F value Pr(>F)
Species 2 69.47 34.74 26.07 5.41e-10 ***
Residuals 110 146.57 1.33
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# muscle N by species
muscleN.aov <- aov(Muscle.d15N ~ Species, data = mydata)
summary(muscleN.aov)
Df Sum Sq Mean Sq F value Pr(>F)
Species 2 210.8 105.38 56.49 <2e-16 ***
Residuals 110 205.2 1.87
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# muscle C/N ratio by species
muscleCN.aov <- aov(log(Muscle.CN) ~ Species, data = mydata)
summary(muscleCN.aov)
Df Sum Sq Mean Sq F value Pr(>F)
Species 2 0.437 0.2184 2.044 0.134
Residuals 110 11.757 0.1069
Identifying which species are different from each other can be acheived using post-hoc tests that correct the p-values to take account of the multiple testing involved. Just looking at the d13C muscle data in this example (you can try the others as you like):
TukeyHSD(muscleC.aov)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Muscle.d13C ~ Species, data = mydata)
$Species
diff lwr upr p adj
Pacific.cod-Arrowtooth.Flounder 2.0063584 1.3338843 2.6788325 0.0000000
Pollock-Arrowtooth.Flounder 0.8457546 0.1926797 1.4988296 0.0073897
Pollock-Pacific.cod -1.1606038 -1.7606068 -0.5606007 0.0000343
The larger Pacific Cod appears to be the largest of the fish on average, and they also have the larget d13C and d15N values. We might rightly want to know if this is due to their larger body size (length) or whether it is true species variation. We use analysis of covariance to investigate.
# specify a color pallete with 3 colours for this example
palette(c("black", "red", "blue"))
par(mfrow = c(1,2))
plot(Muscle.d13C ~ Length, col = mydata$Species, data = mydata,
main = "Muscle d13C")
legend("topleft", levels(mydata$Species), col = 1:3, pch = 1)
plot(Muscle.d15N~ Length, col = mydata$Species, data = mydata,
main = "Muscle d15N")
legend("topleft", levels(mydata$Species), col = 1:3, pch = 1)
# reset the palette
palette("default")
d13C.length.model <- glm(Muscle.d13C ~ Species + Length, data = mydata)
summary(d13C.length.model)
Call:
glm(formula = Muscle.d13C ~ Species + Length, data = mydata)
Deviance Residuals:
Min 1Q Median 3Q Max
-3.8335 -0.3881 0.3021 0.7082 1.7101
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -2.024e+01 2.867e-01 -70.593 < 2e-16 ***
SpeciesPacific.cod 1.927e+00 2.981e-01 6.464 2.96e-09 ***
SpeciesPollock 7.635e-01 2.913e-01 2.621 0.010 *
Length 4.975e-04 5.781e-04 0.861 0.391
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for gaussian family taken to be 1.335639)
Null deviance: 216.05 on 112 degrees of freedom
Residual deviance: 145.58 on 109 degrees of freedom
AIC: 359.31
Number of Fisher Scoring iterations: 2
d15N.length.model <- glm(Muscle.d15N ~ Species + Length, data = mydata)
summary(d15N.length.model)
Call:
glm(formula = Muscle.d15N ~ Species + Length, data = mydata)
Deviance Residuals:
Min 1Q Median 3Q Max
-2.7003 -1.0161 0.1646 0.8171 3.4259
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.375e+01 3.152e-01 43.629 < 2e-16 ***
SpeciesPacific.cod 2.910e+00 3.277e-01 8.880 1.52e-14 ***
SpeciesPollock 5.111e-01 3.202e-01 1.596 0.113
Length 2.707e-03 6.354e-04 4.260 4.35e-05 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for gaussian family taken to be 1.613847)
Null deviance: 415.96 on 112 degrees of freedom
Residual deviance: 175.91 on 109 degrees of freedom
AIC: 380.69
Number of Fisher Scoring iterations: 2
While d15N appears to scale positively with length (and you can probably think why this might be), d13C does not appear to be affected by length. You could assess this by comparing a model without length using AIC, where we are looking for AIC to be lower by more than 2 units if we are to justify the inclusion of length in the model:
d13C.species.model <- glm(Muscle.d13C ~ Species, data = mydata)
summary(d13C.species.model)
Call:
glm(formula = Muscle.d13C ~ Species, data = mydata)
Deviance Residuals:
Min 1Q Median 3Q Max
-3.6913 -0.4828 0.2836 0.7510 1.6583
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -20.0782 0.2144 -93.668 < 2e-16 ***
SpeciesPacific.cod 2.0064 0.2830 7.088 1.37e-10 ***
SpeciesPollock 0.8458 0.2749 3.077 0.00264 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for gaussian family taken to be 1.332492)
Null deviance: 216.05 on 112 degrees of freedom
Residual deviance: 146.57 on 110 degrees of freedom
AIC: 358.08
Number of Fisher Scoring iterations: 2
Based on this comparison, with an AIC of 359.31 for the model with length, and an AIC of 358.08 without length, we can conclude that length does not add sufficient explanatory information to warrant its inclusion.
Task: You might want to check whether there is an interaction between species and length in this relationship. Essentially this involves fitting 3 lines, one for each species by length. You can add an interaction to the glm() formulae, using a : to specify which variables interact.
We might well ask whether the isotope values are the same between the liver and muscle tissues within the individual fish. To acheive this, we can conduct a pair-wise t-test that adds power over a standard t-test by treating the observations as being matched, or paired. Since we can’t use the formula method (y~x) easily when we have two columns (vectors) of data to compare, we can’t use the data = mydata option, and so we have to use the $ notation to provide the two vectors for comparison. We can’t even use the subset = option so instead we have to subset each of our vectors manually: mydata$Muscle.d13C[mydata$Species == "Arrowtooth.Flounder"].
You could take this code, copy it and modify it to run paired t-tests on the other species.
t.test(mydata$Muscle.d13C[mydata$Species == "Arrowtooth.Flounder"],
mydata$Liver.d13C.crtd[mydata$Species == "Arrowtooth.Flounder"],
paired = TRUE)
Paired t-test
data: mydata$Muscle.d13C[mydata$Species == "Arrowtooth.Flounder"] and mydata$Liver.d13C.crtd[mydata$Species == "Arrowtooth.Flounder"]
t = -6.1193, df = 28, p-value = 1.328e-06
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-1.8971660 -0.9455698
sample estimates:
mean of the differences
-1.421368
Visualising the multivariate data is easiest if we add ellipses over each species group as we did in the first practical session. The code that follows is taken from there and modified accordingly to make sure the \(x\) and \(y\) data names are appropriate for this new dataset.
Now we can make a nice plot with ggplot taking care to change the aesthetics to make our column names if we just copied and pasted this in from the first practical file.
# load the library
library(ggplot2)
# use our ellipse function to generate the ellipses for plotting
# this is the basic plot, which im not actually going to plot
# as i won't call print(first.plot). Doing this is purely choice,
# and for me it keeps my code tidier and easier to interpret
first.plot <- ggplot(data = mydata, aes(Muscle.d13C, Muscle.d15N)) +
geom_point(aes(color = Species), size = 2)+
ylab(expression(paste(delta^{15}, "N (\u2030)")))+
xlab(expression(paste(delta^{13}, "C (\u2030)"))) +
theme(text = element_text(size=15))
# decide how big an ellipse you want to draw
# NB 50% ellipses this time for no reason other than i dont need them huge
# to get a sense for their size and shape.. indeed we could plot
# Standard Ellipses using p.ell <- stats::pchisq(1, df = 2) which
# results in 0.39
p.ell <- 0.50
# create our plot based on first.plot above
# adding the stat_ellipse() geometry. We
# specify thee ellipse to be plotted using
# the polygon geom, with fill and edge colour
# defined by Taxon as a grouping variable,
# using the normal distribution and with
# a quite high level of transparency.
ellipse.plot <- first.plot +
stat_ellipse(aes(group = Species,
fill = Species,
color = Species),
alpha = 0.3,
level = p.ell,
type = "norm",
geom = "polygon")
print(ellipse.plot)
Analysing these data to test whether the means of both the d13C and the d15N muscle data are simultaneously different can be achieved using Multivariate Analysis of Variance (MANOVA) or PERMANOVA if you are not happy with the assumptions of MANOVA, e.g. that the data are multivariate normal and that they have similar (assumed the same) covariance. In this case, both assumptions seem pretty reasonable to me so I will go with MANOVA here in the first instance. Both models require that we bind the two (or more) column vectors of data together into a matrix which we can do using cbind(). Remember that both these methods take the subet = option should you wish to restrict your comparison to a subset of groups. A useful trick here might be to use the match function (which has an odd, but nice to use alias %in%) and in this example something like: subset = (Species %in% c("Arrowtooth.Flounder", "Pacific.cod")) could provide a useful template for you, especially should you wish to restrict a larger dataset to only the fish or only the invertebrates if they were all in the same dataset. Alternateively, you could create new data.frames of only the data you want, and call them something helpful; you could acheive this with the function dplyr::filter()
multivar.model <- manova(cbind(Muscle.d13C, Muscle.d15N) ~ Species,
data = mydata)
summary(multivar.model)
Df Pillai approx F num Df den Df Pr(>F)
Species 2 0.53826 20.253 4 220 3.223e-14 ***
Residuals 110
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Running the PERMANOVA analysis requires the package vegan so you will need to install.packages("vegan") if you do not already have it. The PERMANOVA method is a randomisation process that jumbles up the data among the groups and determines how likely our observed data are given random chance. It therefore does make the same parametric assumptions about the data being multivariate normal distributed with a common covariance structure among groups.
library(vegan)
Loading required package: permute
Loading required package: lattice
This is vegan 2.5-7
# extract the isotope data for muscle for
# subsequent modelling as a response variable.
Y_muscle <- with(mydata, cbind(Muscle.d13C, Muscle.d15N))
# run a PERMANOVA model
perm.model <- adonis(Y_muscle ~ Species,
data = mydata,
method = "euclidean",
permutations = 9999)
# print output of the permanova model to screen
perm.model
Call:
adonis(formula = Y_muscle ~ Species, data = mydata, permutations = 9999, method = "euclidean")
Permutation: free
Number of permutations: 9999
Terms added sequentially (first to last)
Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
Species 2 280.23 140.117 43.815 0.4434 1e-04 ***
Residuals 110 351.77 3.198 0.5566
Total 112 632.01 1.0000
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Some additional context we might like to add to this analysis is how far apart the centroids are. The centroids of these ellipses are simply the mean(d13C) and the mean(d15N) of each group, and we have already calculated these for the graph above and the results held in the object summarise.muscle. We can then extract the two columns of means, and use the dist() function in R to calculate the pairwise euclidean distances between them.
centroids <- cbind(summarise.muscle$musC, summarise.muscle$musN)
dist.btw.centoids <- dist(centroids)
print(dist.btw.centoids)
1 2
2 3.898797
3 1.278276 2.651899
# calculate the euclidean distances between observations
dis <- vegdist(Y_muscle,
method ="euclidean")
# analysis of multivariate homogeneity of
# group dispersions (variances)
mod <- betadisper(dis, mydata$Species)
# print the results of the ANOVA to screen
knitr::kable(anova(mod))
| Df | Sum Sq | Mean Sq | F value | Pr(>F) | |
|---|---|---|---|---|---|
| Groups | 2 | 7.336073 | 3.6680365 | 4.79097 | 0.0101156 |
| Residuals | 110 | 84.217605 | 0.7656146 | NA | NA |
# perform Tukey's posthoc test between the means
mod.HSD <- TukeyHSD(mod)
# plot the post hoc test results
par(mfrow = c(1,1))
plot(mod.HSD)
As per the help files for ?capscale and the related ?dbrda, if we use Euclidean distances, as we do here, then capscale is identical to rda with the latter also being more efficient. Depending on the details of what one wants to achieve with this analysis which has its origins in PRIMER, you may want to explore more about these alternative functions.
# fit rda model from the vegan package
CAP_like_analysis <- vegan::rda( Y_muscle ~ Species,
data = mydata,
dist = "euclidean")
# try to predict a value for an "Arrowtooth.Flounder"
# for example.
predict(CAP_like_analysis,
newdata = data.frame(Species =
"Arrowtooth.Flounder"))
Muscle.d13C Muscle.d15N
1 -20.07823 14.64197